Thermomechanical properties of a single hexagonal boron nitride sheet m o (N Oh! < (N (N a o > ^' o m X Sandeep Kumar Singh^, M. Neek-Amal^^^*, S. Costamagna^"'^, F. M. Peeters ^ ^Department of Physics, University of Antwerp, Groenenborgerlaan 171, B-2020 Antwerpen, Belgium '^Department of Physics, Shahid Rajaee University, Lavizan, Tehran 16785-136, Iran ^Facultad de Ciencias Exactas Ingenieria y Agrimensura, Universidad Nacional de Rosario and Instituto de Fisica Rosario, Bv. 27 de Fehrero 210 his, 2000 Rosario, Argentina. (Dated: April 23, 2013) Using atomistic simulations we investigate the thermodynamical properties of a single atomic layer of hexagonal boron nitride (h-BN). The thermal induced ripples, heat capacity, and thermal lattice expansion of large scale h-BN sheets are determined and compared to those found for graphene (GE) for temperatures up to 1000 K. By analyzing the mean square height fluctuations {h'^) and the height-height correlation function H{q) we found that the h-BN sheet is a less stiff material as compared to graphene. The bending rigidity of h-BN: i) is about 16% smaller than the one of GE at room temperature (300 K), and ii) increases with temperature as in GE. The difference in stiffness between h-BN and GE results in unequal responses to external uniaxial and shear stress and different buckling transitions. In contrast to a GE sheet, the buckling transition of a h-BN sheet depends strongly on the direction of the applied compression. The molar heat capacity, thermal expansion coefficient and the Gruneisen parameter are estimated to be 25.2 JmoP^K"'^, 7.2xlO~®K~^ and 0.89, respectively. I. INTRODUCTION A single layer of hexagonal boron-nitride (h-BN) is a wide gap insulator that is very promising for opto- electronic technologies [l|, 0], tunnel devices and field- effect transistors [SHa]. According to the well known Mermin- Wagner theorem Q, the stability of any two dimensional crystal is only possible in the presence of atomic corrugations which distort the perfect honey- comb lattice and allow the atoms to explore the out-of- plane direction. Experimental observations have found ripples in suspended sheets of GE [3, Q and atomistic simulations suggest that the strong bonds between the carbon atoms determine the features of these ripples [9j . Understanding the behavior of the ripples is important for many reasons [Ifll- They affect the electronic trans- port properties, e. g., in GE the ripples are believed to be one of the dominant scattering sources which limits the electron mobility [ll|, [l2] ■ h-BN ribbons doped by carbon has recently been pro- posed [3, [l^. In addition, BN based nanostructures are potential materials for thermal management applica- tions |14l4l3 | because of their high thermal conductivity and sensitivity to isotopic substitution and etc. There- fore, the knowledge of the shape and the temperature dependence of the intrinsic ripples is fundamental to de- vise novel nano-devices based on this material. Both GE and h-BN sheets have a honeycomb lat- tice structure, however the different electronic proper- ties and phonon band structure |20l - |22| results in un- equal morphologies and corrugations. Transmission elec- tron microscopy is widely used to resolve the individ- * corresponding author: ncckanial@srttu.edu Armchair O Boron O Nitrogen FIG. 1: (Color online) Schematic view of the single h-BN sheet. Smaller-yellow (bigger-blue) circles refer to the B (N) atoms. The rectangles indicate the atoms that are fixed dur- ing compression. Dashed (straight) lines correspond to arm- chair (zig-zag) uniaxial compression in the direction given by the arrows. Open arrows indicate the shear stress applied in the armchair direction. ual atoms in suspended h-BN sheets [23| where ripples inherently exist. First-principle calculations have been performed using small unit cells, periodically replicated, which are unable to model long wavelength corrugations which require thousands of atoms [2J| while the mechan- ical properties of a h-BN sheet can be estimated by us- ing DFT 123]. Conversely, atomistic simulations using molecular dynamics simulations (MD) enable to study thermo-mechanical properties directly on large scale sam- ples. The modified TersofF potential [26,] (TP) (parame- terized originally for carbon and silicon) with various set of parameters have shown to predict reasonable values for the thermo- mechanical properties of the h-BN sheet. In the pioneer work by Albe et al re-parameterized TP was used to study the impact of energetic boron and nitrogen atoms on a cubic-BN target 27]. Some other groups have also parameterized TP using various experimental data, e.g. Sekkal et al [2a] treated h-BN as a one-component system, using the same potential parameters for both boron (B) and nitrogen (N) (neglecting the B-N interac- tion) to investigate the structural properties. Matsunaga et al Q proposed the TP of B in order to simulate cu- bic boron carbonitrides which are cornpatible with the parameters of nitrogen fitted by Kroll [30| , and recently, Liao et al [31| and Sevik et al [32| reported TP param- eters that were used to study the mechanical properties and the thermal conductivity of a h-BN sheet, respec- tively. In this study, we investigate the thermal rippling be- havior of free standing monolayer h-BN by using state of the art molecular dynamics (MD) simulations. We adopted the TP potential re-parametrized by Sevik et al which has been shown to represent the experimen- tal structure and the phonon dispersion of the two- dimensional h-BN sheet. We found that h-BN is more corrugated and a less stiff material as compared to GE. The height-height correlations can be explained by the theory for continuum membranes [33]. In addition, we report results of both uniaxial and shear stress of a h-BN sheet and compare it with those found for GE. The buck- ling transition for compressed h-BN occurs earlier than for GE and the induced pattern of ripples when subjected to either uniaxial or shear stress shows significant differ- ences. This paper is organized as follows. In Sec. II, we in- troduce the atomistic model and the simulation method. Then, in Sec. Ill we analyze the behavior of the thermal ripples of a h-BN sheet. Here, we obtain the bending rigidity, the heat capacity, the thermal expansion coeffi- cients, and we study the buckling transition of the h-BN sheet when uniaxial external strain and shear stress are applied. All the results are compared with the ones ob- tained for a GE sheet. Finally, in Sec. IV, we present our conclusions. II. METHODS Classical atomistic molecular dynamics (MD) simula- tion is employed to simulate large scale samples of h-BN and GE at temperatures varying from 10 K up to 1000 K. We used a modified Tersoff potential (which is an ordi- nary defined potential in the LAMMPS package [33.l35|) according to the parameters proposed by Sevik et al for the h-BN sheet. All the parameters and a detailed de- scription of the potential energy is given in Ref . [32| . To simulate GE we have used the AIREBO potential [3g which is suitable for simulating hydrocarbons. We em- ployed NPT ensemble simulations with P=0 using the Nose-Hoover thermostat which enables us to allow the size of the system to equilibrate (i.e. the size of the system is not fixed). All the reported values have been computed averaging over 300 — 400 snapshots taken over uncorrelated samples. We start with a square shaped h-BN sheet with periodic boundary conditions and initial dimensions 322A x32ll (315A xSlsA for GE) in the x and y direc- tion, respectively, which correspond to a total number of A^=37888 atoms, and which are sufficiently large in order to capture the long-wavelength regime. Periodic boundary conditions were adopted in both directions. To analyze the thermal ripples we computed the diffraction pattern which is typically studied in experi- ments to detect the size and shape of the corrugations [3| . We obtained also the mean square value of the out-of plane displacements JJh^) of the atoms and, by follow- ing previous works [33, |3g], the height- height correla- tion function {H(q)) which was determined by consid- ering an average of the height between the first neigh- bors of each atom. The latter was shown to follow a q^^ behavior that is expected from the theory of con- tinuum two-dimensional membranes at large values of q in the harmonic approximation (see below). To analyze the differences between the strain induced corrugations in the h-BN and the GE sheets we applied uniaxial and shear stress on both systems as is schematically shown in Fig. [T] In order to apply strain we set the force on the two ends equal to zero and move the end atoms with an infinitesimal compression step {dx = O.OlA) in the desired direction. After each compression step we wait for 2 ps to allow the system to relax and to ensure that the system is in thermal equilibrium [39|. Uniaxial com- pressive stress is applied along the zig-zag or armchair direction separately, and the shear stress is applied on the armchair edges with the opposite velocity for the top and bottom edges. The details of the used method of ap- plying the boundary stress can be found in our previous studies [H-Iill The TP function [H] used in the LAMMPS pack- age [SJ, I33 can be expressed as E e^.4e ^^j) with i) = Yl Yl Mrrj)[fRinj) + bijfAirrj)] (1) (2) j(>«) where /c, /r and /a are cutoff functions, the repulsive pair term, and the attractive pair term, respectively. Ty and bij are respectively the distance from atom i to atom j and the bond order function. The use of TP disregard contributions coming from charge re-distribution which may occur in an ionic crystal. The inclusion of this effect in h-BN modifies the phonon spectrum significantly only for energies corresponding to the optical modes [42|, |43| . The large scale thermal ripples addressed here are gov- erned mainly by the transversal acoustic mode (TA), (c) 0.1 n 1 1 i — -3 -2 -1 0.05 -0.05 -0.1 -0.1 -0.05 0.05 0.1 FIG. 2: (Color online) (a) Contour plot of the heights for an arbitrary snapshot taken during the MD simulation at 300 K. (b) Corresponding simulated diffraction pattern for the sample shown in (a), (c) Zoom view of (b) around q=(0,0). ,\\ ' I ' I ' I (a) GE h-BN 200 400 600 800 1000 T [K] 0.125 0.25 0.5 q [A"'] FIG. 3: (Color onhne) Variation of (a) < h^ > and (b) the bending rigidity versus temperature for both h-BN (open squares) and GE (solid circles), (c) Height-height correlation function H{q) of h-BN at two different temperatures as is indicated. The dashed lines are the harmonic q~* law. which accounts for out-of plane fluctuations, and it cou- ples with the in-plane modes which renormalizes the long wave-length behavior, e.g. the bendiiig rigidity k can be calculated directly from the TA mode l4J] . Therefore, the charge redistribution is expected not to affect the ther- mal fluctuations analyzed here and the use of the TP is justified ^^. III. RESULTS AND DISCUSSION A. Thermal excited ripples In Fig.[31Ja) we show a height contour plot of the atoms of the h-BN sheet for an arbitrary snapshot taken during the MD simulation at 300 K. The corresponding modeled diffraction pattern is shown in Fig.[2][b). This pattern is very similar to the one obtained for the GE sheet Q with the main difference in the distance between the Bragg points due to the unequal lattice constant of h-BN and GE. From the zoom plot around q=(0,0) (Fig.[2][c)) one observes the local structure of the central Bragg point for these intrinsic thermal ripples of h-BN. Notice that the lack of the presence of the q=(0,0) component is a consequence of the absence of a perfectly flat h-BN sheet. The signatures of the thermal induced ripples can also be seen in the mean square value of the out-of-plane fluctuations {h"^)- In Fig. ^a.) we show (/i^) as func- tion of temperature. In comparison with GE (included also in this figure for comparative purposes) we observe that (/i^) is always larger for the h-BN sheet. The weaker (stronger) the atomic B-N (C-C) bonds, the larger (smaller) the corrugations in h-BN (GE). Notice that the B-N bond is not pure covalent and it has an ionic charac- ter which is due to the different electronegativity between the two elements, i.e. the electrons are localized closer to the N atoms rather than the B atoms. Although this partially ionic structure of the h-BN layer increases the interlayer interaction resulting in a larger hardness of 3D bulk h-BN relative to graphite, it makes the single layer of BN less stiff than graphene. Moreover, the stacking of h-BN sheets in bulk h-BN is AAA stacking which is different for the ABC stacking in graphite j23i] . B. Bending rigidity K Accordingly to the two-dimensional theory of contin- uum membranes the height-height correlation function, in the harmonic approximation, where the out-of-plane and in-plane modes are decoupled, is expected to behave as iJ(g) = (|%)|2} = NkuT kSoq A ' (3) where k is the bending rigidity of the membrane, iV is the number of atoms of the sample, 5*0 the surface area per atom and fcs is the Boltzmann constant. In the large wavelength limit, i.e. for g — > 0, the height fluctuations are suppressed by anharmonic couplings between bend- ing and stretching modes giving rise to a renormalized q- dependent bending rigidity and hence Eq. ([3]) is no longer valid [Hill 0- In Fig. ^h) we show k for h-BN cal- culated from the harmonic part of H(q) between q=0.5 A~^ and q=l A~^. In agreement with the larger values of (/i^), we observe that the h-BN membrane posses a smaller k as compared to GE and in the whole temper- ature range it is about 16 % smaller at room tempera- ture (300 K). Note that k for both h-BN and GE increase with temperature. In Fig. |3Kc) we show H{q) at 200 K and 1000 K were the harmonic behavior can be clearly observed (fitted with a dashed line) and, as expected, with increasing temperature H{q) is shifted to larger q. This figure also reveals that the ripples are not charac- terized by an unique wave-length and instead follows the behavior expected from continuum membrane theory. Before ending this section it is worthwhile to investigate an alternative method to estimate the bending rigidity (flexural rigidity). A common method for calculating the bending rigidity of BN-sheet is by performing simulations of BN- nanotubes as a function of its radius (R) of the curved tubes and then extrapolating the results to R^> oo. Hence, one can calculate the elastic en- ergy of the nanotube as a function of the inverse square of the radius, E ~ ^kR^^. This method was used in Ref. [25] and in our previous work [41] to calculate the zero temperature bending rigidity of GE and h-BN which were found to be 1.46 eV and 1.29 eV, respectively. The result obtained with the Tersoff potential using Eq. (3) is less than the result of Ref. [25]. The bending rigidity of GE is larger than h-BN with about 0.15 eV in agree- ment with Ref. [25]. In order to have an indepen- dent check we estimated the bending rigidity of a BN sheet by performing several ground state simulations for (n,n) BN-nanotubes with n=5,.., 20 using the Tersoff potential. Figure 4 shows the variation of the strain energy per atom as func- tion of the inverse square radius of the tube. Fit- ting a line to the data and dividing by the area of half of unit cell atom results in k = 0.86 eV. The result of Tersoff potential agrees well with our 0.00 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 FIG. 4: (Color online) Variation of strain energy versus inverse square of BN-nanotube radius using Tersoff potential. 10 200 400 600 800 1000 T[K] FIG. 5: (Color online) Variation of (a) the total energy per atom and (b) the averaged bond length versus temperature for h-BN and GE (The error bar is less than 10~* eV/atom for the total energy). estimation for k using Eq. (3) (Ref. [25]). Thus we can conclude that result based on the Tersoff potential underestimate the bending rigidity as compared to the result from DFT. C. Heat capacity, thermal expansion and Gruneisen parameter The variation of the total energy per atom with tem- perature (> lOK) is shown in Fig. [5ja). The total en- ergy varies linearly with temperature and gives the corre- sponding lattice contribution to the molar heat capacity at constant pressure (the average size of the system after relaxation is taken constant) which for the h-BN sheet resuhs into Cp = dEr ~dT = 25.2(3) J mor^iC-^ (4) which is comparable to the one for a GE sheet, i.e. 24.5(9) J mol~^K~^ and close to the well known classical molar heat capacity, i.e., C = 3K ~ 24.94 J mol~^K^^, i.e., the Dulong-Petit value, where SR is the universal gas constant. Notice that the heat capacity is a tem- perature dependent parameter that will decrease for temperature below the Debye temperature. However, the classical MD simulation gives the correct high temperature asymptotic limit i.e., the Dulong-Petit value, but fails in the low and intermediate temperature range. In Fig. [5jb) we show the variation of the averaged B- N bond length with temperature. The linear behavior enables us to calculate the linear thermal expansion of the lattice parameter of a h-BN sheet as 1 da a dT lBN = -^ = 7.2(1) xlO-^K-' (5) where a= 1.442 A is the zero temperature lattice parame- ter. The obtained 7 is 33% larger than the one measured for cubic BN, i.e. 4.8xlO~^K~^ [43| and is comparable to the one found for GE, i.e. jge = 10.0(7) x IQ-^ K-\ The latter (i.e. "(nE) is in good agreement with our pre- vious studies |39l.l40|. From the obtained values of 7 and C we can compute the Gruneisen parameter asN = -pr ^ 0.89, Cp (6) where B is the two dimensional bulk modulus and p is the mass density. Using Bh^BN=SeVk^'^ [Bqe = 12.7 eVA-2 HI)), ph-BN = 2A.82/Sh-BN - 3.81 x lO^-^gm-^ [pQj^ = 7.6 X 10-4gm-2) for h-BN (GE), we obtain ubn = 0.89 (aGB=1.2). Note that the estimated value of aoE is better than the one found in Ref. |49|, i.e. -0.2, and is closer to the experimental result, i.e. 2.0 |50| . D. Buckling transition The different stiffness between the h-BN and GE sheets results in different buckling transitions. To study this, we compressed the systems uniaxially where we considered both zig-zag and armchair directions. This was done by fixing one row of atoms in each edge during the compres- sion steps as is indicated schematically in Fig. [TJ The compression rate was taken p = 0.5m/ s which is simi- lar to the one used in our previous study [40|, and small enough to guarantee that the system is in equilibrium during the whole compression process |5l| . The simula- tions were carried out at room temperature. Figure [BJ^a) shows the variation oi < h^ > versus strain, i.e. e = pt/l 0.01 0.1 strain (%) (c) 150 100 50 < -50 llHWUn r'l i 1 -100 -150 ;,< ,t f ,t t ,f 1 7 6 5 4 3[A] 2 1 -1 -2 -3 -4 -100 0. x(A) 100 FIG. 6: (Color online) Variation of < ft^ > versus external uniaxial strain in the (a) armchair and (b) zig-zag direction. Contour plot of the heights of compressed h-BN sheet samples subjected to a compressive strain for fixed < h^ >= 20 A^ for compression in (c) zig-zag direction of the h-BN sheet. Arrows indicate the stress direction. where t is the time (after starting the compression) and I is the initial box length in the direction of the compres- sion. The buckling transition for the h-BN sheet kept at room temperature is found to be 0.6 (0.1) which is smaller than the one for GE, i.e. 1.0 (1.2), for uniaxial compression along the zig-zag (armchair) direction 52]. Hence, when the compression is applied in the zig-zag direction, the h-BN sheet resists much more against the external uniaxial stress as compared to the case the stress is applied along the armchair direction. The smaller critical strain at which the buckling transition in the armchair direction of the h-BN sheet takes place indicates its more sensitive nature to exter- nal uniaxial stress [53]. Although the same argument holds for GE, the difference between the two directions is much smaller. Notice that DFT calculations re- sult in a deviation in bending rigidity (flexural rigidity) between ZZ BN-nanotube and AC BN- 0.65 0.0001 , ,0.0002 U' (A)"' 0.0003 FIG. 7: (Color online) Variation of buckling strain with the longitudinal length of the BN sheet which is compressed along the zig-zag direction. nanotube (ZZ becomes larger than AC) for radius smaller than ~ SA [25] while for larger radius they are the same. A contour plot of the buckled h-BN sheet with (ft,^)=20 A^ and compressed in the zig-zag direction is shown in Fig. IHc). The obtained buckling transition for GE, i.e. 0.8%, agrees very well with our previous studies [39| and is in the range of the experimental measured buckling transition for suspended GE, i.e. 0.5%-0.7% ^E^. It is also interesting to compare the buckling strain with those predicted by the theory of loaded beam. Euler theory for a two end-loaded beam having length L pre- dicts that efc oc L~^ [40, lla^. The first demonstration by MD of Euler buckling in nanostructures goes back to the pioneer work by Yakobson et al |.57J . We performed several simulations in order to find the length dependence of the buckling strain for BN sheets which are e.g. compressed along the zig-zag direction. Figure [7] shows the variation of ei, with the longitudinal length (L) of the zig-zag BN sheet. The solid line is a fit to L~^ and the symbols are the results of our MD sim- ulations which are qualitatively in agreement with the theory of loaded beam. Finally, we report the results for h-BN and GE sheets under shear stress. Here we applied the stress only in the armchair direction as is described schematically in Fig.[T] We found that, due to the larger stiffness, to reach the same value of {h"^) a larger shear stress has to be applied in GE as compared to h-BN. In Figs. |8l^a,b) one can ob- serve significant differences between the corrugated con- figurations of h-BN and GE sheets. These samples were subjected to a shear stress of e = 1.5%. While h-BN is highly sensitive to shear stress and deforms easily, GE can resist larger stress values due to its larger rigidity. IV. CONCLUSIONS The thermal properties of a boron nitride sheet were studied using large scale atomistic simulations. We showed that the scaling properties of a h-BN sheet fol- lows closely the results of membrane theory and hence the thermal excited ripples are not characterized by any particular wave-length. Using the harmonic part of the height-height correlation function we found an increas- ing bending rigidity with temperature which is smaller than the one of graphene. We found that the buckling transition for h-BN depends on the applied compression direction and is much smaller than the one of graphene. The obtained molar heat capacity agrees very well with the well-known Dulong-Petit number, 25.2 J mol~^K^^ and the thermal expansion coefficient was found to be positive and equal to 7.2 xlQ~^K~^. The Gruneisen parameter 0.89 is found to be smaller than the one for graphene, i.e. 1.2. We showed that the different stiffness between the GE and h-BN sheets leads to different pat- terns of deformations in the presence of either uniaxial or shear stress. Acknowledgments We thank K. H. Michel and D. A. Kirilenko for their useful comments on the manuscript. M. N.-A. is supported by the EU-Marie Curie IIF postdoc Fel- lowship/299855. S. Costamagna is supported by the Belgian Science Foundation (BELSPO). This work is supported by the ESF-EuroGRAPHENE project GON- GRAN, the Flemish Science Foundation (FWO-Vl), and the Methusalem programme of the Flemish Government. [1] X. Blase, A. Rubio, S. G. Louie, and M. L. Cohen, Phys. Rev. B 51, 6868 (1995). [2] K. Watanabe, T. Taniguchi, and H. Kanda, Nat. Mater. 3, 404 (2004). [3] L. BritneU, R. V. Gorbachev, R. Jalil, B. D. Belle, F. Schedin, M. I. Katsnelson, L. Eaves, S. V. Morozov, A. S. Mayorov, N. M. R. Peres, A. H. Castro Neto, J. Leist, A. K. Geim, L. A. Ponomarenko, and K. S. Novoselov, Nano Lett. 12, 1707 (2012). [4] G-H. Lee, Y-J. Yu, C. Lee, C. Dean, K. L. Shepard, P. Kim, and J. Hone, Appl. Phys. Lett. 99, 243114 (2011). [5] J. Beheshtian, A. Sadeghi, M. Neek-Amal, K. Michel, and F. M. Peeters, Phys. Rev. B 86, 195433 (2012). [6] N. D. Mermin and H. Wagner, Phys. Rev. Lett. 17, 1133 (1966). [7] J. C. Meyer, A. Geim, M. I. Katsnelson, K. S. Novoselov, 0.8 0.6 0.4 0.2 . [A] -0.2 -0.4 -0.6 -0.8 FIG. 8: Effect of a shear stress of about 1.5 % on the out-of-plane deformation of (a) the h-BN sheet and (b) the graphene sheet. Arrows indicate the stress direction on the armchair edges. Note the larger values of the heights in h-BN. [9 [lo; [11 [12; [is; [14 [is; [16; [it; [is; [19 [2o; [21 [22; [23; [24 T. J. Booth, and S. Roth, Nature (London) 446 60 (2007). D. A. Kirilenko, A. T. Dideykin, and G. Van Tendeloo, Phys. Rev. B 84, 235417 (2011). A. Fasohno, J. H. Los, and M. I. Katsnelson, Nat. Mater. 6, 85 (2007). S. Costamagna, M. Neek-Amal, J. H. Los, and F. M. Peelers, Phys. Rev. B 86, 041408(R) (2012). M. L Katnelson, and A. K. Geim, Phil. Trans. R. Soc. A 366, 195 (2008). S. Costamagna, O. Hernandez, and A. Dobry, Phys. Rev. B 81, 115421 (2010). J. Jung, Z. Qiao, Q. Niu, and A. H. MacDonald, Nano Letters, 12, 2936 (2012). T. Ouyang, Y. Chen, Y. Xie, K. Yang, Z. Bao, and J. Zhong, Nanotechnology 21, 245701 (2010). I. Savic, D. A. Stewart, and N. Mingo, Phys. Rev. B 78, 235434 (2008). D. A. Stewart, L Savic, and N. Mingo, Nano Lett. 9, 81 (2009). C. W. Chang, A. M. Fennimore, A. Afanasiev, D. Okawa, T. Ikuno, H. Garcia, D. Li, A. Majumdar, and A. Zettl, Phys. Rev. Lett. 97, 085901 (2006). H. Shen, Comp. Mater. Sci. 47, 220 (2009). Y. Xiao, X. H. Yan, J. X. Cao, J. W. Ding, Y. L. Mao, and J. Xiang, Phys. Rev. B 69, 205415 (2004). M. Sakurai, Y. Sakai, and S. Saito, J. Phys.: Conf. Ser. 302, 012018 (2011). M. Topsakal, E. Aktiirk, and S. Ciraci, Phys. Rev. B 79, 115442 (2009). D. Golberg, Y. Bando, Y. Huang, T. Terao, M. Mitome, C. Tang, and C. Zhi, ACS Nano 4, 2979 (2010). N. Alem, R. Erni, C. Kisielowski, M. D. Rossell, W. Gan- nett, and A. Zettl, Phys. Rev. B 80, 155425 (2009). H. Sahin, S. Cahangirov, M. Topsakal, E. Bekaroglu, E. Aktiirk, R. T. Senger, and S. Ciraci, Phys. Rev. B 80, 155453 (2009). [25] K. N. Kudin and G. E. Scuseria, and B. L Yakobson, Phys. Rev. B 64, 235406 (2001). [26] J. Tersoff, Phys. Rev. B 37, 6991 (1988). [27] K. Albe and W. MoUer, Comput. Mater. Sci. 10, 111 (1998). [28] W. Sekkal, B. Bouhafs, H. Aourag, and M. Cartier, J. Phys.: Condens. Matter 10, 4975 (1998). [29] K. Matsunaga, C. Fisher, and H. Matsubara, Japan J. Appl. Phys. 39, 48 (2000). [30] P. M. KroU, PhD Thesis (Technische Hochschule Darm- stadt, 1996). [31] M. L. Liao, Y. C. Wang, S. P. Ju, T. W. Lien, and L. F. Huang, J. Appl. Phys. 110, 054310 (2011). [32] C. Sevik, A. Kinaci, J. B. Haskins, and T. Qagin, Phys. Rev. B 84, 085409 (2011). [33] D. Nelson, T. Piran and S. Weinberg, Statistical Mechan- ics of Membranes and Surfaces (Word Scientific, Singa- pore, 2004). [34 [35^ [36 [37 [38; [39 [40 [41 [42; |http://lammps. sandia.gov S. J. Plimpton, J. Comp. Phys. 117, 1 (1995). S. J. Stuart, A. B. Tutein, and J. A. Harrison, J. Chem. Phys. 112, 6472 (2000). S. Costamagna and A. Dobry, Phys. Rev. B 83, 233401 (2011). J. H. Los, M. I. Katsnelson, O. V. Yazyev, K. V. Za- kharchenko, and A. Fasolino, Phys. Rev. B 80, 121405 (2009). M. Neek-Amal and F. M. Peelers, Phys. Rev. B 82, 085432 (2010); Appl. Phys. Lett. 97, 153118 (2010). M. Neek-Amal and F. M. Peelers, J. Phys.: Condens. Matter 23, 045002 (2011). A. Lajevardipour, M. Neek-Amal, and F. M. Peelers, J. Phys.: Condens. Matter 24, 175303 (2012). K. H. Michel and B. Verbeck, Phys. Status Solidi B 246, 2802 (2009). [43] K. H. Michel and B. Verbeck, Phys. Rev. B 83, 115328 (2011). [44] L.J. Karssemeijer and A. Fasolino, Surf. Sci. 605, 1611 (2011). [45] S. S. Han, J. K. Kang, H. M. Lee, A. C. T. van Duin, and W. A. Goddard, J. of Chem. Phys. 123, 114703 (2005). [46] P. Le Doussal and L. Radzihovsky, Phys. Rev. Lett. 69, 1209 (1992). [47] G. Demazeau, Diamond and Related Materials, 2, 197 (1993). [48] K. V. Zakharchenko, M. L Katsnelson, and A. Fasolino, Phys. Rev. Lett. 102, 046808 (2009). [49] V. Perebeinos and J. Tersoff, Phys. Rev. B 79, 241409(R) (2009). [50] T. M. G. Mohiuddin, A. Lombardo, R. R. Nair, A. Bonetti, G. Savini, R. Jalil, N. Bonini, D. M. Basko, C. Galiotis, N. Marzari, K. S. Novoselov, A. K. Geim, and A. C. Ferrari, Phys. Rev. B 79, 205433 (2009). [51] M. Neek-Amal, Phys. Rev. Lett 106, 209701 (2011). [52] A movie of the buckling process in a zig-zag BN sheet is provided in the Supplementary Material. [53] Q. Peng, W. Ji, and S. De, Comp. Mat. Science 56, 11 (2012). [54] O. Frank, G. Tsoukleri, J. Parthenios, K. Papagelis, I. Riaz, R. Jalil, K. S. Novoselov, and C. Galiotis, ACS Nano 4, 3131 (2010). [55] G. Tsoukleri, J. Parthenios, K. Papagelis, R. Jalil, A. C. Ferrari, A. K. Geim , K. S. Novoselov, and C. Galiotis, Small 5, 2397 (2009). [56] J. Singer, J. Arbocz and T. Weller, Experimental Methods in Buckling of Thin- Walled Structures (Wiley, New York, 1998). R. M. Jones, Buckling of Bars, Plates, and Shells (Taylor and Francis, Philadelphia, 1999). [57] B. L Yakobson, C. J. Brabec, and J. Bernholc, Phys. Rev. Lett. 76, 2511 (1996).